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Abstract 

The Standard Vicsek Model (SVM) is a minimal nonequilibrium model of self- 
propelled particles that appears to capture the essential ingredients of critical 
flocking phenomena. In the SVM, particles tend to align with each other and 
form ordered flocks of collective motion; however, perturbations controlled by a 
noise term lead to a noise-driven, continuous order-disorder phase transition. In 
this work, we extend the SVM by introducing a parameter a that allows particles 
to be individualistic instead of gregarious, i.e. to choose a direction of motion 
independently of their neighbors. By focusing on the small-noise regime, we 
show that a relatively small probability of individualistic motion (around 10%) 
is sufficient to drive the system from a Vicsek-like ordered phase to a disor- 
dered phase. Despite the fact that the a— extended Model preserves the 0(n) 
symmetry, the interaction range, as well as the dimensionality of the underlying 
SVM, this novel phase transition is found to be discontinuous (first-order), an 
intriguing manifestation of the richness of the nonequilibrium flocking/ swarming 
phenomenon. 



1. Introduction 

Swarming and flocking arise as emergent phenomena of collective motion be- 
havior in a large variety of living and non-living self-propelled particle systems 
of great interdisciplinary interest. Nature offers a huge variety of collective mo- 
tion phenomena in self-propelled living systems at all scales, from biomolecular 
micromotors, migrating cells and growing bacteria colonies, to insect swarms, 
fish schools, bird flocks, mammal herds and even human crowds Moreover, 
many non-living systems of great practical interest involve collective motion 
and swarming behavior, particularly in robotics, where swarms of robots are 
used in terrain exploration, plague control, optimization of telecommunication 
networks, surveillance and defense, and other tasks without centralized control 
that appear too challenging to be carried out by an individual agent Very 
recently, the novel concept of chemical robots (also known as chobots) has been 
envisioned as an army of millions of micrometer-sized robots whose tasks will be 
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to release a chemical payload, or to mix two chemical reactants from different 
compartments within the chobots when they reach their goal 

Instead of focusing on the specific details that make each of these self- 
propelled systems unique, statistical physicists have been studying the general 
patterns of collective motion, aiming to identify the general laws and underlying 
principles that may govern their behavior. From this perspective, one important 
question to address is the onset of ordered macroscopic phases, i.e. the way in 
which individuals having short-range interactions are capable of self-organizing 
into large-scale cooperative patterns in the absence of leaders or other ordering 
cues from the environment. By analogy with large molecular systems, flocking 
and swarming phenomena can be associated with phase transitions that depend 
on a few parameters that characterize the macroscopic states, such as the den- 
sity of individuals and the flock size. For instance, Buhl et al. [H investigated 
the collective motion of locusts, which display a density-driven transition from 
disordered movement of individuals within the group to highly aligned collec- 
tive motion. Similar transitions have been observed in the colective motion of 
zooplankton swarms fish schools and many other self-propelled particle 
systems ([l[ and references therein). 

On the theoretical side, Vicsek et al. [1] proposed a minimal model to study 
the onset of order in systems of self-driven individuals, which was later fol- 
lowed by other investigations by means of agent- based modeling @, H 0, the 
Newtonian force-equation approach [13, 1 1 1 . and the hydrodynamic approxi- 
mation [IlEIlHl. The so_called Standard Vicsek Model (SVM) [g| assumes 
that neighboring individuals tend to align their direction of movement when 
they are placed within a certain interaction range. This alignment rule, which 
would trivially lead to fully ordered collective motion, is complemented by a 
second one that introduces noise in the communications among individuals. In 
their seminal paper, Vicsek et al. show that the noise amplitude drives the 
system through a continuous, second-order transition between an ordered phase 
of collective motion and a disordered phase. 

In this context, the aim of this work is to explore an extension of the SVM, 
in which an additional parameter a is introduced to control the gregarious 
vs individualistic behavior of particle motion. That is, an individual has a 
probability a to adopt its own direction of motion regardless of its neighbors, 
and a probability 1 — a to move according to the SVM's rules. This a— extended 
Model may account for "free will" behavior in biological systems (i.e. the fact 
that living organisms are not ruled by fixed decision algorithms and are therefore 
able to make unforeseen individualistic decisions at any time) as well as for 
random failures in robotic and other artificial self-propelled particle systems (e.g. 
the chemically-driven chobots mentioned above). By focusing on the small-noise 
regime, we show that a relatively small probability of individualistic motion 
(around 10%) is sufficient to drive the system from a Vicsek-like ordered phase 
to a disordered phase. Besides the practical interest of extending the well- 
studied SVM to novel scenarios of particle behavior, we find a theoretically 
intriguing manifestation of the richness of the swarming phenomenon, namely 
that the a— driven phase transition is discontinuous (first-order), despite the 
fact that the a— extended Model preserves the 0(n) symmetry, the interaction 
range, as well as the dimensionality of the underlying SVM. 

The rest of this paper is organized as follows. In Section 2, we present 
the definition of the SVM and the a— extended Model. Section 3 presents a 
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discussion of our results, while our conclusions finally appear in Section 4. 



2. The Standard Vicsek Model (SVM) and the a-extended Model 

The SVM [6| consists of a fixed number of interacting particles, AT, which are 
moving on a plane. In computer simulations, that plane is typically represented 
by a square of side L with periodic boundary conditions [y, 15]. The particles 
move off-lattice with constant and common speed vq = \v\. Each particle inter- 
acts locally adopting the direction of motion of the subsystem of neighboring 
particles (within an interaction circle of radius Rq centered in the considered 
particle), which is perturbed by the presence of noise. Since the interaction 
radius is the same for all particles, we adopt the interaction radius as the unit 
of length throughout, i.e. Rq = 1. The model is implemented as a cellular 
automaton, i.e. all particles update their states simultaneously in one time 
step. 

The updated direction of motion for the i— th particle, is given by 



gi+l _ 
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vQ , (i) 



where rj is the noise amplitude, the summation is carried over all particles within 
the interaction circle centered at the i— th particle, and £| is a realization of a 
(^-correlated white noise uniformly distributed in the range between — ir and it. 

In this work, we implement the model dynamics by adopting the so-called 
forward update (FU) rule, in which the updated velocities at time t+1 determine 
the new positions according to 

£i +1 =xt + v? +1 M , (2) 

where At = 1 is the unitary time step of the cellular automaton. This scheme 
is used for computational convenience, since the transient length is shorter than 
using the backward update rule that determines the new velocities after the 
positions are updated. In the context of Vicsek model studies, the FU rule was 
first introduced by Chate et al. HI 17 1 . As explained in Ref. , the FU scheme 



requires to implement an algorithm that ensures the O(n) rotational invariancc 
of the SVM, since spurious effects induced by the boundary conditions of the 
simulation space may arise otherwise. At every time step, a frame rotation 
angle, 6/ r , is randomly selected with uniform probability in the — 7r < 0/,- < tt 
range. After all particle velocities are accordingly transformed, the usual Vicsek 
model rules are applied. The effects of boundary conditions on the displacement 
space, as well as alternative methods to ensure the 0(n) rotational invariance, 
are further discussed in Refs. 18, 1^. lipj. 



In order to study the relative effects of gregarious vs individualistic behavior, 
we consider an extension of the SVM model in which a parameter a > is intro- 
duced as follows. When considering the new velocity of a given particle, we first 
determine whether that particle will behave individualistically or gregariously: 
with probability a, the particle will choose a new direction of motion at random 
(and independently of its neighbors), and with probability 1 — a, the particle's 
new velocity will follow the standard Vicsek's rule, Eq.([T]). In the former case, 
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the new direction of motion is randomly chosen with uniform probability in 
the range between —n and tt. In passing, let us point out that, in a different 
context, mixtures of spins coupled to independent heat baths have been studied 
for the kinetic Ising model and far-from-equilibrium spin chains 21 , 2^, 0, [HI ■ 
By analogy, we can consider our model as dealing with off-lattice ferromagnetic 
"spins" that undergo displacements along the direction of the spin. At any 
given time, gregarious individuals could be thought of as "cold spins", while 
individualistic individuals could be analogously regarded as "hot spins" coupled 
to a thermal bath at infinite temperature. Notice also that, although more com- 
plex variations of the SVM could be considered, we purposefully focus on the 
simplest extension of the SVM that includes a "free-will" option (or "failure" 
option, if considered within the context of non-living systems). Indeed, it is 
our aim to preserve the minimal nature of the SVM and show that, from the 
point of view of nonequilibrium statistical mechanics, a qualitatively novel and 
distinct behavior arises from this simple extension. Moreover, the SVM is fully 
recovered by simply letting q->0. 

The parameters of the model are the number of particles N , the dimension 
of the displacement space d, the linear size of the displacement space L, the 
interaction radius Rq, the particle density p, the particle speed vq, the noise 
amplitude rj, and the individualistic probability a. Notice, however, that p = 
N/L d , where d — 2 in the standard case, and that Rq = 1 can be chosen as 
the unit of length. Moreover, v^ 1 merely plays the role of a "thermalization 
parameter" that measures how many times, on average, two neighbors check 
out each other's positions while they remain at a distance within the unitary 
interaction radius. Therefore, relative large values of speed (typically w > 0.3) 
correspond to the low-thermalization regime, characterized by highly anisotropic 
diffusion and the manifestation of simulation artifacts in the form of directionally 
quantized density waves (see Ref. [lj| for details). In order to avoid spurious 
artifacts, we choose v = 0.1 and the density value p = 0.25 throughout, which 
are standard values in the Vicsek model literature @. Moreover, in order to 
pinpoint effects due to the a— extension, we consider the small noise regime 
(rj <C 1) throughout and focus on the behavior of the system over a wide range 
of values of the parameter a. Also, finite-size effects are assessed by studying 
several system sizes in the range 16384 < N < 131072. 

In our simulations, one time step is defined by the simultaneous update of 
velocities and positions for all N particles in the system. For each set of parame- 
ter values, time series of 7 x 10 6 time steps have been generated. All observables 
measured here are statistical averages within the steady-state regime, which is 
typically attained after 1 — 2 x 10 6 time units. Let us also point out that, for the 
largest system sizes studied here (i.e. N > 10 5 ), generating this large number 
of steady-state configurations requires a significant computational effort. 



3. Results and Discussion 

The well-studied SVM exhibits a far-from-equilibrium continuous phase tran- 
sition between ordered states of motion at low noise levels and disordered motion 
at high noise levels. The natural order parameter is given by the absolute value 
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Figure 1: Plots of the order parameter ip versus the parameter a (which measures the prob- 
ability for individualistic behavior) for different system sizes, as indicated. By varying a, 
the system is driven across an order-disorder transition with very mild finite-size effects. The 
dashed straight line shows the mean-field behavior for rj = and p — > oo in the thermodynamic 
limit, where the effective interaction becomes long-range and the order-disorder transition is 
exactly located at at = 1. 



of the normalized mean velocity of the system [6j : 

1 N 



where ip is close to zero in the disordered phase and close to unity in the ordered 
phase. In order to analyze the critical nature of the SVM, this noise-driven phase 
transition has been investigated by several approaches, e.g . finite-size scaling 



analysis and short-time critical dynamic techniques [25 . 

In this work, we explore the a— extended Vicsek model, where the parameter 
a measures the probability for individualistic behavior (see Section 2 above for 
details). Since we focus on the small- noise regime, the a — >• limit corresponds 
to the SVM ordered phase. Figure [1] shows the order parameter <p as a function 
of the parameter a. We observe that a relatively small probability of individ- 
ualistic motion, at ~ 0.09, is sufficient to drive the system from a Vicsek-like 
ordered phase (<p > 0) to a disordered phase (<p « 0). Moreover, finite-size 
effects are mild, as indicated by the near-overlap of the plots of different system 
sizes in the range 16384 < N < 131072. The absence of significant finite-size 
effects near the transition region at at ~ 0.09 suggests that the phase transi- 
tion may be discontinuous (first-order) instead of the continuous (second-order) 
phase transition characteristic of the SVM. Indeed, the Binder cumulant analy- 
sis presented below will provide very strong evidence to confirm the first-order 
nature of the phase transition and will also allow us to determine the transi- 
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Figure 2: Typical snapshot configurations for TV = 16384 and different values for a, namely: 

(a) a = 0.02, (b) a = 0.093, and (c) a = 0.20. These configurations have been generated 
by applying the algorithm proposed in Ref. Ilgj . in which the 0(n) rotational invariance is 
enforced explicitly in the model's dynamic rules. 

tion value with great accuracy. The dashed straight line in Figure Q] shows the 
mean-field limit of (p vs a for r\ — when all particles interact with each other, 
i.e. when p — > oo in the thermodynamic limit (iV — > oo). In this limiting case, 
a fraction a of the system is formed by particles moving in random directions 
and therefore has a null contribution to the order parameter, while a fraction 
1 — a is fully aligned along the mean direction of motion of the whole system. 
Thus, cp = 1 — a. It is interesting to note that, in this limit of infinite density 
where the effective interaction becomes long-range, the transition from order to 
disorder is exactly located at at = 1. 

Further qualitative insight can be gained by looking at the clustering of par- 
ticles in typical snapshot configurations across the transition. Notice that the 
standard order parameter for Vicsek-like flocking phenomena, Eq.([3]), refers to 
the overall direction of motion of the particles. Therefore, the spatial distribu- 
tion of particles offers a related, complementary view of the degree of order in the 
system. Figure [2] shows snapshot configurations for a system of size N — 16384 
and different values for the individualistic probability a, namely: (a) a = 0.02, 

(b) a — 0.093, and (c) a = 0.20. The low-a regime (Fig. [5Ja)) is characterized 
by the existence of large flocks that are responsible for a system- wide, large-scale 
coherent motion. Inbetween flocks, some regions display relatively large voids in 
which the particle density is very low. At the transition (Fig.[2jb)), the ordered 
flocks appear scattered throughout the system in coexistence with the disor- 
dered phase. Notice also that flock sizes do not appear to depart significantly 
from the average flock size, i.e. the snapshot suggests that the flock size distri- 
bution does not follow a power law (as would be expected if the transition were 
continuous). Finally, the high-a regime (Fig. [2jc)) shows a configuration where 
individuals are nearly evenly distributed across the system and large flocks are 
absent, as expected for the disordered phase. Thus, at a qualitative level, we 
observe that the spatial distribution of individuals exhibits different configura- 
tions that, as the parameter a is increased, change from ordered to disordered 
through a state of phase coexistence. Quantitative evidence for this first-order 
phase transition is provided by the so-called Binder cumulant of Vicsek's order 
parameter, as follows. 
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Figure 3: Binder fourth-order cumulant U4 as a function of the individualistic probability a 
for different system sizes, as indicated. 



The Binder cumulant, denned by 



'-'-If' 



is a fourth-order cumulant dependent on the variance and the kurtosis of the 
order parameter probability distribution. Since, for second-order phase transi- 
tions, the scaling prefactor of the cumulant is independent of the sample size, 
plots of U4 versus the control parameter lead to a common (size-independent) 
intersection point that corresponds to the location of the critical value of the or- 
der parameter in the thermodynamic limit [27| . Moreover, for first-order phase 
transitions, a characteristic signature is for U4 to have a sharp fall towards neg- 
ative values 0, [2^]. This hallmark behavior is somewhat elusive because it 
only appears for relatively large system sizes. 

Figure [3] displays the Binder cumulant as a function of the individualis- 
tic probability a for different system sizes, as indicated. In agreement with 
previous figures, the cumulant's behavior shows the existence of an a— driven 
phase transition. Furthermore, although not apparent in the small system size 
(N = 16384), the results for the larger systems studied here (N > 65536) show 
sharp minima at the transition point, which, as discussed above, is the signa- 
ture behavior of discontinuous (first-order) phase transitions. The sharpness of 
these minima allows us to determine accurately the location of the transition at 
a t = 0.0925 ±0.010. 

A hallmark behavior of a discontinuous phase transition is the phenomenon 
of metastability. In equilibrium systems, this phenomenon results from the 
coexistence of two or more local minima of the free energy, where each minimum 
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Figure 4: Time series of the order parameter ip for a system of size N = 131072 and different 
values of a across the transition region, (a) For a = 0.086, the system is ordered and the 
order parameter displays small fluctuations around (p 0.25 — 0.30. (b) For a = 0.0928, 
strong excursions are observed between one ordered (meta)stable state at ip fa 0.20 and one 
disordered (meta)stable state at ip fa 0.05. (c) For a = 0.093, the ordered metastable state 
is less likely than the disordered one. (d) For a = 0.0947, the system stays chiefly in the 
disordered state; excursions to the ordered metastable state become shorter in duration and 
less frequent. 

corresponds to one of the coexisting phases at the transition. In finite-size 
systems, fluctuations may drive the system from one local minimum to another 
one across the free energy barrier separating those (meta)stable states. In our 
nonequilibrium model, an analogous behavior is observed in time series of the 
order parameter. Figure U shows the order parameter ip for a system of size 
N = 131072 as a function of time for different values of a across the transition 
region. Fig. [Ha ) corresponds to a = 0.086, i.e. below the phase transition. The 
system is pinned at the ordered phase and displays relatively small fluctuations 
around ip « 0.25 — 0.30. Notice also the transient regime from an initially 
fully disordered state a t = to the steady-state regime, which in this case is 
reached after i«5x 10 5 . Fig. 0|b) corresponds to a — 0.0928, where strong 
excursions are observed between one ordered (meta)stable state at ip w 0.20 
and one disordered (meta)stable state at ip w 0.05. By increasing slightly the 
control parameter to a — 0.093, Fig. @Jc) shows that the ordered metastable 
state becomes less likely than the disordered one. This strong sensitivity to 
small variations in a is in agreement with the sharp transition displayed by the 
Binder cumulant (Fig. [3] above). By increasing the control parameter further to 
a = 0.0947, as shown by Fig. EJd), the system stays mostly in the disordered 
state. In fact, excursions to the ordered metastable state become noticeably 
shorter in duration and much more sporadic. 

The phenomenon of metastability is also reflected in the behavior of the prob- 
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Figure 5: Probability distribution functions (pdf) of the order parameter for a system of size 
N = 131072 and different values of a across the transition region, (a) For a = 0.086, the 
system is ordered and the pdf is unimodal and peaked around tp 0.25 — 0.30. (b) For 
a = 0.0928, the distribution is bimodal and corresponds to the coexistence of two distinct 
phases, (c) For a = 0.093, the peak corresponding to the disordered phase prevails, although 
the distribution is still bimodal. (d) For a = 0.0947, the peak associated with the ordered 
phase has vanished and the distribution is completely dominated by the peak associated with 
the disordered phase. 

ability distribution functions (pdf) of the order parameter across the transition 
region. Figure [5] shows the order parameter probability distribution functions 
for a system of size N = 131072 and the same values of a used in the previous 
Figure. In Fig. [SJ^a), the distribution for a — 0.086 is unimodal and centered 
around ip ~ 0.25 — 0.30, i.e. the order parameter has small fluctuations around 
the positive mean value that characterizes the ordered phase. Fig. [5jb) , which 
corresponds to the transition value a = 0.0928, is a bimodal distribution with 
two similar peaks associated with the coexistence of the ordered and disordered 
phases. In Fig. [Sfc), which was obtained with a = 0.093, one observes that the 
peak corresponding to the disordered phase becomes dominant, although the 
distribution is still bimodal. Finally, Fig. [51(d) represents the disordered regime 
for a = 0.0947, where the local maximum associated with the ordered phase 
has completely vanished. As a summary, the panels of Figure [5] display the 
occurrence of a sharp order-disorder transition that takes place within a very 
narrow range of values of a. The hallmark of a discontinuous phase transition, 
in which the two phases coexist, becomes evident by the bimodal probability 
distribution function of the order parameter at the transition, i.e. a « ott- Be- 
low the transition within the ordered phase (a < a t ), the peak corresponding 
to the disordered phase vanishes; and viceversa, above the transition within the 
disordered phase (a > at), the peak associated with the ordered phase vanishes 
and the peak corresponding to the disordered phase prevails. 
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Figure 6: Binder fourth-order cumulant U4 as a function of the individualistic probability a 
for N = 131072 and different noise amplitudes, as indicated. Arrows mark the occurrence 
of minima. Inset: phase diagram in the r\ vs a parameter space, where filled circles indicate 
first-order transitions, the star shows the pure-Vicsek limit with a second-order transition, 
and the red diamond shows a hypotetical tricritical point. Sec more details in the text. 



So far, we have focused on the effects of adding individualistic motion, as 
controlled by the parameter a, on Vicsek-like self-propelled particles within 
the small- noise regime (rj <C 1). Let us now investigate the behavior of the 
a— extended Model under non-negligible noise amplitudes. Figure [5] shows the 
Binder cumulant as a function of the individualistic probability a for N — 
131072 and different noise amplitudes, as indicated. For very small noise am- 
plitudes, e.g. 77 = 0.01, the Binder cumulant shows a sharp minimum at the 
transition, which as we have discussed above, is a characteristic signature of 
first-order phase transitions (recall Fig. [3]). As the noise amplitude is increased, 
however, the minima become less noticeable and appear shifted towards smaller 
values of a. Increasing the noise amplitude beyond the pure-Vicsek critical 
transition, i.e. rj > rj c = 0.134 [25| . all signatures of a transition vanish and 
the Binder cumulant appears fiat, no longer depending on the parameter a. 
Based on these findings, we are able to sketch a phase diagram in the 77 vs a 
parameter space, as displayed in the Inset to Fig. [5] On the one hand, transition 
points with characteristic first-order signatures are shown by filled circles. On 
the other hand, the collective motion transition in the pure-Vicsek limit with 
a = has been characterized as second-order |(|; for the velocity and density 
parameters used in this work, the continuous second-order phase transition has 
been located at r\ c = 0.134 [25| (shown by a star in the Inset to Fig. [S]). 

Certainly, the fact that the smooth a— extension of the SVM changes the 
critical nature of the system is very intriguing. Indeed, these findings are a 
display of the richness of the swarming phenomenon even when represented by 
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means of a minimal nonequilibrium model such as Vicsek's. As an example of a 
rich equilibrium minimal model, in which the critical nature of the underlying 
system is modified by the addition of a Hamiltonian term that preserves the 
model symmetries, the interaction range and the dimensionality, let us mention 
the Blume-Capel extension of the classical Ising model [30l . l31| . The Blumc- 
Capel model adds a symmetry-preserving term proportional to a crystal field 
D that favors impurities (S = states) over up/down spins (S — ±1 states). 
In the D — s- — oo limit, the Blume-Capel model tends smoothly to the classical 
two-state Ising model. However, in some regions of the T — D parameter space, 
the order-disorder phase transition becomes discontinuous (first-order) and a 
tricritical point can be determined. Reasoning by analogy, and based on the 
results discussed above, we therefore conjecture the existence of a tricritical 
point (as indicated in the Inset to Fig. [5]), which would separate the standard 
noise-driven continuous transition from the a— driven discontinuous transition 
reported here. The existence of such phenomenon, however, remains an open 
question that lies beyond the scope of the present work. 

4. Conclusions 

In this work, we investigated the so-called a— extended Vicsek Model, in 
which a parameter a was introduced to control the gregarious vs individualistic 
behavior of particle motion. By focusing on the small-noise regime, we showed 
that a relatively small probability of independent particle motion was sufficient 
to drive the system from a Vicsek-like ordered phase to a disordered phase. 
Indeed, the order parameter and its associated fourth-order Binder cumulant 
showed very strong evidence of a discontinuous, first-order phase transition oc- 
curring at at — 0.0925 ± 0.010. Further qualitative and quantitative insight 
into the nature of the transition was provided by snapshot configurations, time 
series of the order parameter, as well as probability distribution functions of the 
order parameter. 

The aims of the present study were two-fold. On the one hand, by extending 
the well-studied Standard Vicsek Model (SVM) to novel scenarios of particle be- 
havior (while preserving, at the same time, its nature as a minimal model), we 
may account for important features of real self-propelled particle systems that 
are not captured by the SVM, as for instance free-will behavior in biological sys- 
tems and random failures in robotic and other artificial systems. Moreover, on 
the other hand, we found a theoretically intriguing manifestation of the richness 
of the swarming phenomenon, namely that the a— driven phase transition was 
discontinuous (first-order), despite the fact that the a— extended Model pre- 
served the 0(n) symmetry, the interaction range, as well as the dimensionality 
of the underlying SVM. 

By analogy to equilibrium systems such as the Blume-Capel extension of 
the classical Ising model, we conjecture the existence of a tricritical point in 
the rj — a parameter space of the a— extended Model, which would separate the 
standard noise-driven continuous transition from the a— driven discontinuous 
transition reported here. Indeed, the existence of such phenomenon remains a 
very interesting open question that certainly deserves further investigation. 
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